Symcap
Teegan Innis
December 29, 2016
Project Summary
This project aimed to look at the spatial variability of color morph and Symbiodinium clades C and D present in the Kane’ohe Bay, O’ahu, Hawai’i population of Montipora capitata. We investigated the distributions at scales ranging from location within the bay to location on an individual reef. We also looked at differences among reef types (fringing vs. patch) and depths. Heterogeneous mixtures of symbiont clades were considered in the analysis for spatial patterns. By investigating spatial variability of Symbiodinium, we furthered the understanding of potential stress-response in Kane’ohe Bay.
Library Packages
knitr::opts_knit$set(root.dir = normalizePath(".."))
setwd("~/Symcap")
library(data.table)
library(lsmeans)
library(devtools)
library(plyr)
library(reshape2)
library(popbio)
library(RgoogleMaps)
library(plotrix)
library(zoo)
library(rgdal)
library(car)
library(scales)
library(png)
library(pixmap)
library(ecodist)
library(cluster)
library(fpc)
library(clustsig)
library(mapplots)
library(foreign)
library(nnet)
library(ggplot2)
library(mlogit)Import/Merge Collection and qPCR Data
Coral_Data <- read.csv("Data/Collection Data/Coral_Collection.csv")
Coral_Data$Depth..m. <- as.numeric(as.character(Coral_Data$Depth..m.))
source_url("https://raw.githubusercontent.com/jrcunning/steponeR/master/steponeR.R")
Mcap.plates <- list.files(path = "Data/qPCR_data", pattern = "txt$", full.names = T)
Mcap <- steponeR(files = Mcap.plates, delim="\t",
target.ratios=c("C.D"),
fluor.norm = list(C=2.26827, D=0),
copy.number=list(C=33, D=3))
Mcap <- Mcap$result
Mcap <- Mcap[grep("+", Mcap$Sample.Name, fixed=T, invert = T), ]
Mcap <- Mcap[grep("NTC", Mcap$Sample.Name, fixed = T, invert = T), ]
Mcap <- Mcap[grep("PCT", Mcap$Sample.Name, fixed = T, invert = T), ]
colnames(Mcap)[which(colnames(Mcap)=="Sample.Name")] <- "Colony"
Mcap$fail <- ifelse(Mcap$C.reps < 2 & Mcap$D.reps < 2, TRUE, FALSE)
fails <- Mcap[Mcap$fail==TRUE, ]
Mcap <- Mcap[which(Mcap$fail==FALSE),]
Mcap$C.D[which(Mcap$C.reps<2)] <- -Inf
Mcap$C.D[which(Mcap$D.reps<2)] <- Inf
Mcap <- Mcap[with(Mcap, order(Colony)), ]
Mcap$propC <- Mcap$C.D / (Mcap$C.D+1)
Mcap$propD <- 1-Mcap$propC
Mcap$propD[which(Mcap$C.D==-Inf)] <-1
Mcap$propC[which(Mcap$C.D==-Inf)] <-0
Mcap$propD[which(Mcap$C.D==Inf)] <-0
Mcap$propC[which(Mcap$C.D==Inf)] <-1
Mcap$Dom <- ifelse(Mcap$propC>Mcap$propD, "C", "D")
Symcap<-merge(Coral_Data, Mcap, by="Colony", all=T)
Symcap <- Symcap[with(Symcap, order(Colony)), ]
Symcap$Mix <- factor(ifelse(Symcap$propC>Symcap$propD, ifelse(Symcap$propD!=0, "CD", "C"), ifelse(Symcap$propD>Symcap$propC, ifelse(Symcap$propC!=0, "DC", "D"), NA)), levels = c("C", "CD", "DC", "D"))
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")Change Bay Area for 18, 26 and F7-18
Symcap$Bay.Area[which(Symcap$Reef.ID=="26")] <- "Central"
Symcap$Bay.Area[which(Symcap$Reef.ID=="18")] <- "Southern"
Symcap$Bay.Area[which(Symcap$Reef.ID=="F7-18")] <- "Southern"Adjust Depth by Mean Sea Level
To account for the influence of tides, depth was adjusted according to the differences in mean sea level from the recorded sea level on each collection day to the nearest 6-minute interval. Mean sea level was obtained from NOAA daily tide tables for Moku o Lo’e.
JuneTide=read.csv("Tide Tables/Station_1612480_tide_ht_20160601-20160630.csv")
JulyTide=read.csv("Tide Tables/Station_1612480_tide_ht_20160701-20160731.csv")
AugustTide=read.csv("Tide Tables/Station_1612480_tide_ht_20160801-20160812.csv")
Tide<-rbind(JuneTide, JulyTide, AugustTide)
Tide$Time <- as.POSIXct(Tide$TimeUTC, format="%Y-%m-%d %H:%M:%S", tz="UTC")
attributes(Tide$Time)$tzone <- "Pacific/Honolulu"
Symcap$Time2 <- as.POSIXct(paste(as.character(Symcap$Date), as.character(Symcap$Time)),
format="%m/%d/%y %H:%M", tz="Pacific/Honolulu")
Symcap$Time=Symcap$Time2
# Add estimated times for missing values
Symcap$Time[170] <- as.POSIXct("2016-06-14 12:07:00")
Symcap$Time[177] <- as.POSIXct("2016-06-14 12:20:00")
Symcap$Time[178] <- as.POSIXct("2016-06-14 12:22:00")
Symcap$Time[180] <- as.POSIXct("2016-06-14 13:08:00")
Symcap$Time[187] <- as.POSIXct("2016-06-14 12:42:00")
Symcap$Time[188] <- as.POSIXct("2016-06-14 12:44:00")
Symcap$Time[206] <- as.POSIXct("2016-06-16 13:10:00")
Symcap$Time[208] <- as.POSIXct("2016-06-16 13:24:00")
Symcap$Time[211] <- as.POSIXct("2016-06-16 12:37:00")
Symcap$Time[218] <- as.POSIXct("2016-06-16 12:27:00")
Symcap$Time[448] <- as.POSIXct("2016-07-16 13:32:00")
Round6 <- function (times) {
x <- as.POSIXlt(times)
mins <- x$min
mult <- mins %/% 6
remain <- mins %% 6
if(remain > 3L) {
mult <- mult + 1
} else {
x$min <- 6 * mult
}
x <- as.POSIXct(x)
x
}
Symcap$Time.r <- Round6(Symcap$Time)
Tide$Time.r <- Tide$Time
merged<-merge(Symcap, Tide, by="Time.r", all.x=T)
merged$newDepth <- merged$Depth..m.- merged$TideHTCoral Collection Map
Collection sites represented on this map indicate the 9 fringe and 16 patch reef sites for a total of 25 collection reefs. In total, 707 colonies were collected for the analysis.
KB <- c(21.46087401, -157.809907)
KBMap <- GetMap(center = KB, zoom = 13, maptype = "satellite", SCALE = 2, GRAYSCALE = FALSE)
save(KBMap, file = "KBMap.Rdata")
load("KBMap.Rdata")
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
XY <- subset(XY, Reef.ID!="37")
par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude, col=153, pch=21, bg="#FF7F50", lwd=2, cex = 1.2)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)
par(new=T, mar=c(14,21,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.87, 21.41, -157.75, 21.52)
box()
Symbiont Dominance and Color Morph Across Geographic Scales
Significant Differences
Dominant Symbiont per Reef Area
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Dom, Symcap$Reef.Area)
chisq.test(results)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 136.26, df = 1, p-value < 2.2e-16prop.table(results, margin = 2)##
## Slope Top
## C 0.7767857 0.3294574
## D 0.2232143 0.6705426par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Reef Area", ylab = "Symbmiont Proportion")
legend("topright", legend=c("C", "D"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)
Dominant Symbiont per Color Morph
results=table(Symcap$Dom, Symcap$Color.Morph)
chisq.test(results)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 284.23, df = 1, p-value < 2.2e-16prop.table(results, margin = 2)##
## Brown Orange
## C 0.95327103 0.32987013
## D 0.04672897 0.67012987par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Color Morph", ylab = "Symbiont Proportion")
legend("topright", legend=c("C", "D"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)
Symbiont Community Composition per Color Morph
results=table(Symcap$Mix, Symcap$Color.Morph)
chisq.test(results)##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 287.29, df = 3, p-value < 2.2e-16prop.table(results, margin = 2)##
## Brown Orange
## C 0.819314642 0.290909091
## CD 0.133956386 0.038961039
## DC 0.043613707 0.649350649
## D 0.003115265 0.020779221par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray50", "gray85", "gray100"), xlab = "Color Morph", ylab = "Symbiont Community Composition")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray50", "gray85", "gray100"), inset = c(-.2, 0), xpd = NA)
Symbiont Community Composition per Reef Area
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Mix, Symcap$Reef.Area)
chisq.test(results)##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 138.97, df = 3, p-value < 2.2e-16prop.table(results, margin = 2)##
## Slope Top
## C 0.678571429 0.275193798
## CD 0.098214286 0.054263566
## DC 0.214285714 0.651162791
## D 0.008928571 0.019379845par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray50", "gray85", "gray100"), xlab = "Reef Area", ylab = "Symbiont Community Composition")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray50", "gray85", "gray100"), inset = c(-.2, 0), xpd = NA)
Symbiont Community Composition per Dominant Symbiont
results=table(Symcap$Mix, Symcap$Dom)
chisq.test(results)##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 707, df = 3, p-value < 2.2e-16prop.table(results, margin = 2)##
## C D
## C 0.86635945 0.00000000
## CD 0.13364055 0.00000000
## DC 0.00000000 0.96703297
## D 0.00000000 0.03296703par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray 10", "gray 85", "gray 40", "gray100"), xlab = "Dominant Symbiont", ylab = "Symbiont Mixture Proportion")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray85", "gray40", "gray100"), inset = c(-.2, 0), xpd = NA)
When looking at D-only colonies, the Ct values are on the higher end (35 or greater) indicating poor amplification or the potential for C amplification being pushed back in the cycle. This is supported by the fact that 5 of the 9 D-only colonies had C present in 1 qPCR replicate.
Proportion of D When Present in Mixture
propD <- merged$propD[which(merged$propD > 0 & merged$propD < 1)]
hist(propD, xlab = "Proportion of Clade D", ylab = "Number of Samples", main = "", col = "gray75")
Color Morph per Reef Area
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Color.Morph, Symcap$Reef.Area)
chisq.test(results)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 127.59, df = 1, p-value < 2.2e-16prop.table(results, margin = 2)##
## Slope Top
## Brown 0.6169265 0.1744186
## Orange 0.3830735 0.8255814par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Reef Area", ylab = "Color Morph Proportion")
legend("topright", legend=c("Brown", "Orange"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)
Color Morph per Bay Area
results=table(merged$Color.Morph, merged$Bay.Area)
prop.table(results, margin = 2)##
## Central Northern Southern
## Brown 0.5096154 0.4153846 0.4440789
## Orange 0.4903846 0.5846154 0.5559211chisq.test(results)##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 3.8811, df = 2, p-value = 0.1436Color Morph per Reef ID
results=table(Symcap$Color.Morph, Symcap$Reef.ID)
chisq.test(results)##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 45.314, df = 24, p-value = 0.005347Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propCol=table(Symcap$Color.Morph, Symcap$Reef.ID)
propCol=prop.table(propCol, margin = 2)
propCol <- as.data.frame.matrix(propCol)
props <- data.frame(t(propCol))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
rownames(XY) <- XY$Reef.ID
XY <- XY[, -1]
XY <- na.omit(XY)
apply(XY, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["Orange"], reef["Brown"]), radius = 7, col = c("orange", "sienna"))
})
Mantel Test for Spatial Autocorrelation
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propCol=table(Symcap$Color.Morph, Symcap$Reef.ID)
propCol=prop.table(propCol, margin = 2)
propCol <- as.data.frame.matrix(propCol)
props <- data.frame(t(propCol))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
reef.dists <- dist(cbind(XY$Longitude, XY$Latitude))
col.dists <- dist(XY$Brown)
set.seed(12456)
mantel(col.dists~reef.dists)## mantelr pval1 pval2 pval3 llim.2.5% ulim.97.5%
## -0.08204377 0.83500000 0.16600000 0.30800000 -0.12641256 -0.02236124Color Morph per Bay Area Adjusted for Depth
To discount the vertical spatial plane (depth), these results are adjusted to “assume” that all colonies were collected from the same depth, therefore considering simply the coordinate plane as the distributional influence on the resulting values.
merged$Color <- ifelse(merged$Color.Morph=="Orange", 1, 0)
dat <- aggregate(data.frame(prop=merged$Color), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T)
mod <- glm(Color ~ newDepth + Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "Raw", "O.Adj")
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propCol=table(Symcap$Color.Morph, Symcap$Reef.ID)
propCol=prop.table(propCol, margin = 2)
propCol <- as.data.frame.matrix(propCol)
props <- data.frame(t(propCol))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
XY2 <- merge(XY, res, by = "Reef.ID")
XY2$B.Adj <- 1-XY2$O.Adj
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["O.Adj"], reef["B.Adj"]), radius = 7,
col = c("orange", "sienna"))
})
Mantel Test for Spatial Autocorrelation
reef.dists <- dist(cbind(XY$Longitude, XY$Latitude))
col.dists <- dist(XY2$O.Adj)
set.seed(12456)
mantel(col.dists~reef.dists)## mantelr pval1 pval2 pval3 llim.2.5% ulim.97.5%
## -0.07099365 0.79400000 0.20700000 0.38600000 -0.11365492 -0.02840648Symbiont Community Composition per Bay Area
results=table(Symcap$Mix, Symcap$Bay.Area)
prop.table(results, margin = 2)##
## Central Northern Southern
## C 0.50961538 0.61538462 0.49174917
## CD 0.08173077 0.05641026 0.09900990
## DC 0.39423077 0.30256410 0.40594059
## D 0.01442308 0.02564103 0.00330033chisq.test(results)##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 14.719, df = 6, p-value = 0.02256Non-Significant Differences
Dominant Symbiont per Reef Type
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 0.77593, df = 1, p-value = 0.3784Dominant Symbiont per Bay Area
results=table(Symcap$Dom, Symcap$Bay.Area)
chisq.test(results)##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 3.8852, df = 2, p-value = 0.1433Dominant Symbiont per Reef ID
results=table(Symcap$Dom, Symcap$Reef.ID)
chisq.test(results)##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 34.598, df = 24, p-value = 0.07459Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propDom=table(Symcap$Dom, Symcap$Reef.ID)
propDom=prop.table(propDom, margin = 2)
propDom <- as.data.frame.matrix(propDom)
props <- data.frame(t(propDom))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
rownames(XY) <- XY$Reef.ID
XY <- XY[, -1]
XY <- na.omit(XY)
apply(XY, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["C"], reef["D"]), radius = 7, col = c("blue", "red"))
})
Mantel Test for Spatial Autocorrelation
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propDom=table(Symcap$Dom, Symcap$Reef.ID)
propDom=prop.table(propDom, margin = 2)
propDom <- as.data.frame.matrix(propDom)
props <- data.frame(t(propDom))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
reef.dists <- dist(cbind(XY$Longitude, XY$Latitude))
dom.dists <- dist(XY$C)
set.seed(12456)
mantel(dom.dists~reef.dists)## mantelr pval1 pval2 pval3 llim.2.5% ulim.97.5%
## 0.15974663 0.04000000 0.96100000 0.04500000 0.08875833 0.24328773pval1 indicates that there is a positive correlation between distance and difference among reefs in terms of dominant symbiont clade. Reefs that are closer in proximity to each other are more similar in terms of symbiont dominance. A positive mantelr value indicates a positive correlation.
Clustering of Similar Reefs
Though insignificant differences resulted when investigating symbiont dominance across reefs, reefs closer in proximity appear more similar than reefs further apart. Clustering was tested here.
doms <- hclust(dom.dists)
plot(doms, labels = XY$Reef.ID)
a <- hclust(reef.dists*dom.dists)
plot(a, labels = XY$Reef.ID)
Dominant Symbiont Plots
Raw Values
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
XY <- subset(XY, Reef.ID!="37")
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")
XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["raw.c"], reef["raw.d"]), radius = 7,
col = c("red", "blue"))
})
Depth Adjustment
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")
XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["c.adj"], reef["d.adj"]), radius = 7,
col = c("red", "blue"))
})
Depth Interaction Adjustment
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")
XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["c.adj.int"], reef["d.adj.int"]), radius = 7,
col = c("red", "blue"))
})
Color Morph and Dominant Symbiont Interaction Plots
A multinomial logistic regression was performed on the interaction of color morph and dominant symbiont to discount the influence of depth on spatial distribution.
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID),
FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")
XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c
manteltable = table(merged$Dom, merged$Color.Morph, merged$Reef.ID)
nc <- aggregate(interaction(merged$Color.Morph, merged$Dom),
by=list(merged$Reef.ID), FUN=table)
nc <- data.frame(Reef.ID=as.character(nc[,1]), prop.table(nc[,2], margin=1))
XY3 <- merge(XY2, nc, by="Reef.ID", all=T)
merged$ColDom <- interaction(merged$Color.Morph, merged$Dom)
model <- multinom(ColDom~newDepth+Reef.ID, merged)
means <- lsmeans(model, specs = "Reef.ID")
means
pp <- fitted(model)
newdat <- data.frame(Reef.ID = levels(merged$Reef.ID),
newDepth = mean(merged$newDepth, na.rm=T))
pred <- predict(model, newdata = newdat, "probs")
new <- data.frame(Reef.ID=as.character(newdat[,1]), pred)
XY4 <- merge(XY3, new, by="Reef.ID", all=T)
PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)
XY4 <- merge(XY3, new, by="Reef.ID", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Reef.ID
XY4 <- XY4[, -1]
XY4 <- na.omit(XY4)
apply(XY4, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["Brown.C.x"], reef["Orange.C.x"], reef["Brown.D.x"],
reef["Orange.D.x"]), radius = 7,
col = c("turquoise4", "steelblue1", "tomato", "coral"))
})
legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("turquoise4", "steelblue1", "tomato", "coral"))
par(new=T, mar=c(15,22,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.9, 21.41, -157.75, 21.53)
box()
Raw proportion values for Color Morph and Dominant Symbiont at each reef not discounting the effect of depth.
Mantel Test for Spatial Autocorrelation
reef.dists <- dist(cbind(XY4$Longitude, XY4$Latitude))
dom.dists <- bcdist(cbind(XY4$Brown.C.x, XY4$Orange.C.x, XY4$Brown.D.x, XY4$Orange.D.x))
set.seed(12456)
mantel(dom.dists~reef.dists)## mantelr pval1 pval2 pval3 llim.2.5% ulim.97.5%
## 0.12673388 0.07400000 0.92700000 0.12300000 0.03226547 0.19014964PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)
XY4 <- merge(XY3, new, by="Reef.ID", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Reef.ID
XY4 <- XY4[, -1]
XY4 <- na.omit(XY4)
apply(XY4, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["Brown.C.y"], reef["Orange.C.y"], reef["Brown.D.y"],
reef["Orange.D.y"]), radius = 7,
col = c("turquoise4", "steelblue1", "tomato", "coral"))
})
legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("turquoise4", "steelblue1", "tomato", "coral"))
par(new=T, mar=c(15,22,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.9, 21.41, -157.75, 21.53)
box()
Adjusted values of Dominant Symbiont and Color Morph adjusted for the influence of depth.
Mantel Test for Spatial Autocorrelation
reef.dists <- dist(cbind(XY4$Longitude, XY4$Latitude))
dom.dists <- bcdist(cbind(XY4$Brown.C.y, XY4$Orange.C.y, XY4$Brown.D.y, XY4$Orange.D.y))
set.seed(12456)
mantel(dom.dists~reef.dists)## mantelr pval1 pval2 pval3 llim.2.5% ulim.97.5%
## 0.05339855 0.22800000 0.77300000 0.49600000 -0.02379201 0.10887683Color and Dominant Symbiont per Bay Area
When considering the effect of bay area on the interaction of color morph and dominant symbiont, the proportion of Brown colonies dominated by D increases as you move from the north to the south of the bay. The proportion increases almost 6x, yet the small number of colonies makes this non-compelling.
Type=table(merged$ColDom, merged$Bay.Area)
prop.table(Type, margin = 2)##
## Central Northern Southern
## Brown.C 0.48076923 0.41538462 0.41254125
## Orange.C 0.11057692 0.25641026 0.17821782
## Brown.D 0.02884615 0.00000000 0.02970297
## Orange.D 0.37980769 0.32820513 0.37953795chisq.test(Type)## Warning in chisq.test(Type): Chi-squared approximation may be incorrect##
## Pearson's Chi-squared test
##
## data: Type
## X-squared = 20.185, df = 6, p-value = 0.002567Color Morph per Reef Type
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: results
## X-squared = 1.5265, df = 1, p-value = 0.2166Symbiont Community Composition per Reef Type
##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 2.8371, df = 3, p-value = 0.4174Symbiont Community Composition per Reef ID
##
## Pearson's Chi-squared test
##
## data: results
## X-squared = 87.127, df = 72, p-value = 0.1081Symbiont Dominance and Color Morph Across Depths
Dominant Symbiont by Depth
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
Dom1 <- subset(merged, !is.na(newDepth) & !is.na(Dominant))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Dominant
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 705 942.15
## newDepth 1 129.01 704 813.13 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1logi.hist.plot(Dom1$newDepth, Dom1$Dominant, boxp = FALSE, type = "hist", col="gray", xlabel = "Depth (m)", ylabel = "", ylabel2 = "")
mtext(side = 4, text = "Frequency", line = 3, cex=1)
mtext(side = 4, text = "C D", line = 2, cex = 0.75)
mtext(side = 2, text = "Probability of clade C Symbiont", line = 3, cex = 1)
merged$DepthInt <- cut(merged$newDepth, breaks = 0:13)
merged$Dominant2 <- ifelse(merged$Dom=="C", 1, 0)
results=table(merged$Dominant2, merged$DepthInt)
results##
## (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
## 0 149 33 27 23 14 4 3 2 2 1 0
## 1 73 52 83 92 46 27 15 12 16 8 5
##
## (11,12] (12,13]
## 0 1 0
## 1 0 1props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(-.2, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Probability of D-Dominance")
Bars indicate relative frequency of occurence of clade C vs. D dominance at 1m depth intervals when pooling all samples collected. Overlaid line indicates probability of a colony being D-dominated across depths.
Mixture by Depth
merged$Mixture <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 1, 0)
merged$Mixture2 <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 0, 1)
results=glm(Mixture~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Mixture
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 705 973.27
## newDepth 1 93.723 704 879.55 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results=table(merged$Mixture2, merged$DepthInt)
results##
## (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
## 0 154 44 39 28 18 8 4 2 5 2 1
## 1 68 41 71 87 42 23 14 12 13 7 4
##
## (11,12] (12,13]
## 0 1 0
## 1 0 1props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Mix", "Non-Mix"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(-.23, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
results=glm(Mixture~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Mixture Proportion")
Bars indicate relative frequency of occurrence of Mixtures vs. Non-Mixtures at 1m depth intervals when pooling all samples collected. The overlaid line indicates the probability of a colony being a mixture across depths.
Color Morph by Depth
merged$Color <- ifelse(merged$Color.Morph=="Orange", 1, 0)
results=glm(Color~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Color
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 706 974.49
## newDepth 1 123.97 705 850.51 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Color <- subset(merged, !is.na(newDepth) & !is.na(Color))
logi.hist.plot(independ = Color$newDepth, depend = Color$Color, type = "hist", boxp = FALSE, ylabel = "", col="gray", ylabel2 = "", xlabel = "Depth (m)")
mtext(side = 4, text = "Frequency", line = 3, cex=1)
mtext(side = 4, text = "Brown Orange", line = 2, cex = 0.75)
mtext(side = 2, text = "Probability of Orange Color Morph", line = 3, cex = 1)
merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results##
## (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
## 0 178 42 48 54 26 8 8 1 2 0 0
## 1 44 43 62 61 34 23 10 13 17 9 5
##
## (11,12] (12,13]
## 0 1 0
## 1 0 1props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.25), alpha("orange", 0.25)), inset = c(-.22, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Probability of Orange-Dominance")
Bars indicate relative frequency of Brown vs. Orange color morph dominance at 1m depth intervals when pooling all samples collected. The overlaid line indicates the probablity of a colony to be Orange across depths.
Color Morph and Dominant Symbiont Interaction
Dominant Symbiont per Color Morph and Reef Area
merged$Reef.Area <- ifelse(merged$Reef.Area!="Top", yes = "Slope", no = "Top")
table(merged$Dom, merged$Color.Morph, merged$Reef.Area)## , , = Slope
##
##
## Brown Orange
## C 263 85
## D 13 87
##
## , , = Top
##
##
## Brown Orange
## C 43 42
## D 2 171model=aov(Dominant~Color.Morph*Reef.Area, data = merged)
Anova(model, type = 2)## Anova Table (Type II tests)
##
## Response: Dominant
## Sum Sq Df F value Pr(>F)
## Color.Morph 41.193 1 317.738 < 2.2e-16 ***
## Reef.Area 5.924 1 45.696 2.907e-11 ***
## Color.Morph:Reef.Area 2.470 1 19.052 1.463e-05 ***
## Residuals 91.011 702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dominant Symbiont per Depth and Bay Area
model=aov(Dominant~newDepth*Bay.Area, data = merged)
Anova(model, type = 2)## Anova Table (Type II tests)
##
## Response: Dominant
## Sum Sq Df F value Pr(>F)
## newDepth 25.026 1 125.1913 < 2e-16 ***
## Bay.Area 0.954 2 2.3870 0.09265 .
## newDepth:Bay.Area 1.555 2 3.8886 0.02092 *
## Residuals 139.933 700
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dominant Symbiont per Color Morph and Depth
Because an interactive effect of reef area and color morph was observed and slope is indicative of a depth gradient, the interaction between depth and color morph was tested here.
model1=lm(Dominant~Color.Morph*newDepth, data = merged)
Anova(model1, type = 2)## Anova Table (Type II tests)
##
## Response: Dominant
## Sum Sq Df F value Pr(>F)
## Color.Morph 46.532 1 368.836 < 2.2e-16 ***
## newDepth 3.496 1 27.708 1.877e-07 ***
## Color.Morph:newDepth 7.347 1 58.238 7.592e-14 ***
## Residuals 88.563 702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Depth at which Orange switches from D to C Dominance
While brown was always C-dominated, orange was observed to switch from D to C dominance. The depth threshold where this switch occurred is represented here.
threshdepth <- function(color) {
df <- subset(merged, Color.Morph==color)
plot(df$Dominant2~df$newDepth, xlab="Depth (m)", ylab = "Probability of Clade C Symbiont",
main=color)
abline(h = 0.5, lty=2)
results=glm(Dominant2~newDepth, family = "binomial", data = df)
pval <- data.frame(coef(summary(results)))$`Pr...z..`[2]
mtext(side=3, text=pval)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted ~ seq(0,11,0.01))
thresh <- ifelse(pval < 0.05,
newdata$newDepth[which(diff(sign(fitted - 0.5))!=0)], NA)
return(thresh)
}
sapply(levels(merged$Color.Morph), FUN=threshdepth)
## Brown Orange
## NA 3.53merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID),
FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")## NOTE: Results may be misleading due to involvement in interactionsres <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")
XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c
manteltable = table(merged$Dom, merged$Color.Morph, merged$Reef.ID)
nc <- aggregate(interaction(merged$Color.Morph, merged$Dom),
by=list(merged$Reef.ID), FUN=table)
nc <- data.frame(Reef.ID=as.character(nc[,1]), prop.table(nc[,2], margin=1))
XY3 <- merge(XY2, nc, by="Reef.ID", all=T)
merged$ColDom <- interaction(merged$Color.Morph, merged$Dom)
model <- multinom(ColDom~newDepth+Reef.ID, merged)## # weights: 108 (78 variable)
## initial value 978.723819
## iter 10 value 654.967406
## iter 20 value 638.634607
## iter 30 value 637.242985
## iter 40 value 636.981882
## iter 50 value 636.866349
## iter 60 value 636.848431
## iter 70 value 636.846645
## final value 636.846634
## convergedmeans <- lsmeans(model, specs = "Reef.ID")
means## Reef.ID prob SE df lower.CL upper.CL
## 10 0.25 NaN 78 NaN NaN
## 13 0.25 NaN 78 NaN NaN
## 14 0.25 1.646361e-10 78 0.25 0.25
## 18 0.25 NaN 78 NaN NaN
## 19 0.25 NaN 78 NaN NaN
## 20 0.25 NaN 78 NaN NaN
## 21 0.25 3.134578e-10 78 0.25 0.25
## 25 0.25 2.842171e-14 78 0.25 0.25
## 26 0.25 2.328306e-10 78 0.25 0.25
## 30 0.25 NaN 78 NaN NaN
## 38 0.25 NaN 78 NaN NaN
## 42 0.25 2.098707e-10 78 0.25 0.25
## 46 0.25 2.851581e-10 78 0.25 0.25
## 5 0.25 NaN 78 NaN NaN
## Deep 0.25 0.000000e+00 78 0.25 0.25
## F1-46 0.25 4.656613e-10 78 0.25 0.25
## F2-25 0.25 NaN 78 NaN NaN
## F3-18 0.25 NaN 78 NaN NaN
## F4-34 0.25 NaN 78 NaN NaN
## F5-34 0.25 NaN 78 NaN NaN
## F6-Lilipuna 0.25 4.115903e-10 78 0.25 0.25
## F7-18 0.25 NaN 78 NaN NaN
## F8-10 0.25 3.292723e-10 78 0.25 0.25
## F9-5 0.25 NaN 78 NaN NaN
## HIMB 0.25 NaN 78 NaN NaN
##
## Results are averaged over the levels of: ColDom
## Confidence level used: 0.95pp <- fitted(model)
depthmod <- multinom(ColDom ~ newDepth, merged)## # weights: 12 (6 variable)
## initial value 978.723819
## iter 10 value 702.413324
## final value 702.105235
## convergeddepthdat <- data.frame(newDepth = seq(0,12,0.5))
depthpred <- predict(depthmod, newdata=depthdat, "probs")
#plot(NA, xlim=c(0,12), ylim=c(0,1))
stackpoly(depthpred[,c(4,2,3,1)], stack=T, col=rev(c("burlywood4","burlywood3","darkorange3","darkorange")),
ylim=c(0,1), xaxlab=depthdat$newDepth, xlab="Depth", ylab="Probability")
legend("top", legend=c("Brown.C", "Brown.D", "Orange.C","Orange.D"), pch=22,
pt.bg=c("burlywood4","burlywood3","darkorange3","darkorange"), inset=-0.3, xpd=NA, bty="n")
df <- subset(merged, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,11,0.01), col="sienna", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Clade D Symbiont", line = 3, cex = 1)
legend("topright", legend=c("Brown", "Orange"), fill=c("sienna", "orange"), inset = c(-.22, 0), xpd = NA)
Dominant Symbiont by Depth and Reef Type
model3=aov(Dominant2~newDepth*Reef.Type, data = merged)
Anova(model3, type = 2)## Anova Table (Type II tests)
##
## Response: Dominant2
## Sum Sq Df F value Pr(>F)
## newDepth 24.785 1 123.5859 < 2.2e-16 ***
## Reef.Type 0.014 1 0.0674 0.795285
## newDepth:Reef.Type 1.643 1 8.1944 0.004327 **
## Residuals 140.785 702
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1df <- subset(merged, Reef.Type=="Patch")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,11,0.01), col="brown1", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Clade D Symbiont", line = 3, cex = 1)
legend("topright", legend=c("Patch", "Fringe"), fill=c("dodgerblue3", "brown1"), inset = c(-.22, 0), xpd = NA)
Color Morph by Depth and Reef Type
model2=aov(Color~newDepth*Reef.Type, data = merged)
Anova(model2, type = 2)## Anova Table (Type II tests)
##
## Response: Color
## Sum Sq Df F value Pr(>F)
## newDepth 27.693 1 133.3174 <2e-16 ***
## Reef.Type 0.074 1 0.3563 0.5508
## newDepth:Reef.Type 1.195 1 5.7553 0.0167 *
## Residuals 146.026 703
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1df <- subset(merged, Reef.Type=="Patch")
results=glm(Color2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
results=glm(Color2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,11,0.01), col="brown1", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Orange Color Morph", line = 3, cex = 1)
legend("topright", legend=c("Patch", "Fringe"), fill=c("dodgerblue3", "brown1"), inset = c(-.22, 0), xpd = NA)
Dominant Symbiont per Color Morph by Depth and Reef Type
Because depth and reef type have an interactive effect on both color morph and dominant symbiont clade, the dominant symbiont per color morph at each reef type was visualized.
df <- subset(merged, Reef.Type=="Patch")
dfo <- subset(df, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = dfo)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Patch")
dfb <- subset(df, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = dfb)
newdata <- list(newDepth=seq(0,11,0.01))
par(new = T, mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="black", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
dfo <- subset(df, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = dfo)
newdata <- list(newDepth=seq(0,11,0.01))
par(new = T, mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="red", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
dfb <- subset(df, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = dfb)
newdata <- list(newDepth=seq(0,11,0.01))
par(new = T, mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
legend("topright", legend=c("PO", "PB", "FO", "FB"), fill=c("dodgerblue3", "black", "red", "orange"), inset = c(-.2, 0), xpd = NA)
mtext(side = 2, text = "Probability of Clade D Symbiont", line = 3, cex = 1)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
Dominant Symbiont by Depth and Reef ID
model2=aov(Dominant~newDepth*Reef.ID, data = merged)
Anova(model2, type = 2)## Anova Table (Type II tests)
##
## Response: Dominant
## Sum Sq Df F value Pr(>F)
## newDepth 23.693 1 122.6080 < 2.2e-16 ***
## Reef.ID 6.905 24 1.4888 0.063088 .
## newDepth:Reef.ID 8.772 24 1.8915 0.006452 **
## Residuals 126.765 656
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1merged$Dominant2 <- ifelse(merged$Dom=="C", 0, 1)
par(mfrow=c(5,5))
domreef <- function(id) {
df <- subset(merged, Reef.ID==id)
results=glm(Dominant2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(1, 1, 3, 1))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l",
col="dodgerblue3", lwd=3, xlab="", ylab="")
mtext(side = 3, text = id)
abline(h=0.5, lty=2)
}
sapply(levels(merged$Reef.ID), FUN=domreef)
On the y-axis, a value of 1 indicates D-dominance and a value of 0 indicates C-dominance.
Mixture by Depth and Reef ID
merged$Mixture <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 1, 0)
model4=aov(Mixture~newDepth*Reef.ID, data = merged)
Anova(model4, type = 2)## Anova Table (Type II tests)
##
## Response: Mixture
## Sum Sq Df F value Pr(>F)
## newDepth 18.913 1 90.9371 < 2.2e-16 ***
## Reef.ID 9.987 24 2.0008 0.003245 **
## newDepth:Reef.ID 8.020 24 1.6068 0.033901 *
## Residuals 136.434 656
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1par(mfrow=c(5,5))
mixreef <- function(id) {
df <- subset(merged, Reef.ID==id)
results=glm(Mixture~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(1, 1, 3, 1))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l",
col="dodgerblue3", lwd=3, xlab="", ylab="")
mtext(side = 3, text = id)
abline(h=0.5, lty=2)
}
sapply(levels(merged$Reef.ID), FUN=mixreef)
On the y-axis, a value of 1 indicates a mixture and a value of 0 indicates single symbiont clade.
Figures
Symbiont Community Composition per Dominant Symbiont
results=table(Symcap$Mix, Symcap$Dom)
chisq.test(results)
prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c(alpha("blue", 0.75), alpha("blue", 0.25), alpha("red", 0.25), alpha("red", 0.75)), xlab = "Dominant Symbiont", ylab = "Symbiont Mixture Proportion")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c(alpha("blue", 0.75), alpha("blue", 0.25), alpha("red", 0.25), alpha("red", 0.75)), inset = c(-.2, 0), xpd = NA)
When looking at D-only colonies, the Ct values are on the higher end of the spectrum (35 or greater) indicating poor amplification and the potential for C amplification being pushed back in the cycle. This is supported by the fact that 5 of the 9 D-only colonies had C present in 1 qPCR replicate.
Proportion of Clade D when Present
propD <- merged$propD[which(merged$propD > 0 & merged$propD < 1)]
propDHist <- subset(merged, propD > 0 & propD < 1)
propDHist$propD <- cut(propDHist$propD, breaks = 5)
DCol=table(propDHist$Color.Morph, propDHist$propD)
z <- with(subset(merged, propD==0), table(Color.Morph))
o <- with(subset(merged, propD==1), table(Color.Morph))
DCol <- cbind(z, DCol, o)
par(mar=c(4, 4, 2, 0))
bars <- barplot(DCol, xlab = "Clade D Proportion", ylab = "Number of Samples", main = "", col = c(alpha("sienna", 0.55), alpha("orange", 0.55)), axisnames = FALSE, space = c(0,0.3,0,0,0,0,0.3))
axis(side = 1, at = c(1.3, 2.3, 3.3, 4.3, 5.3, 6.3), labels = c(">0%", "20%", "40%", "60%", "80%", "<100%"))
axis(side = 1, at = c(0.5, 7.2), labels = c("All C", "All D"), tick = FALSE)
Depth Influence on Dominant Symbiont and Color Morph
par(mfrow=c(3,1))
merged$DepthInt <- cut(merged$newDepth, breaks = 0:13)
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
merged$Dominant2 <- ifelse(merged$Dom=="C", 1, 0)
results=table(merged$Dominant2, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(2, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)),
xlab = "", ylab = "",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(0, 0), xpd = NA)
par(new = T)
par(mar=c(2.1, 4, 2, 6))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="", ylab="Probability of D-Dominance")
merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(3, 4, 1, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)),
xlab = "", ylab = "Probability of Orange-Dominance",
space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.25), alpha("orange", 0.25)), inset = c(0, 0), xpd = NA)
par(new = T)
par(mar=c(3.1, 4, 1, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="", ylab="")
df <- subset(merged, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 0, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="Depth (m)", ylab="Probabilty of D-Dominance", axisnames=FALSE, xaxs = "i", yaxs = "i")
abline(h = 0.5, lty=2)
df <- subset(merged, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0, 11, 0.01), col="sienna", lwd=3)
legend("topright", legend=c("Brown", "Orange"), fill=c("sienna", "orange"), inset = c(0, 0), xpd = NA)
#results=table(merged$Dominant2, merged$Color.Morph)
#chisq.test(results)
#prop.table(results, margin = 2)
#par(new=T, mar=c(10, 10, .5, 6.3))
#barplot(prop.table(results, margin = 2), col = c(alpha("red", 0.35), alpha("blue", 0.35)), xlab = "", ylab = "", yaxt = 'n')
#legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.35), alpha("red", 0.35)), inset = c(0, 0), xpd = NA)Color Morph and Dominant Symbiont at each Reef
par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)
XY4 <- merge(XY3, new, by="Reef.ID", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Reef.ID
XY4 <- XY4[, -1]
XY4 <- na.omit(XY4)
apply(XY4, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["Brown.C.y"], reef["Orange.C.y"], reef["Brown.D.y"],
reef["Orange.D.y"]), radius = 7,
col = c("royalblue3", "cornflowerblue", "red", "rosybrown2"))
})
legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("royalblue3", "cornflowerblue", "red", "rosybrown2"))
par(new=T, mar=c(14,21,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.87, 21.41, -157.75, 21.52)
box()
Adjusted values of Dominant Symbiont and Color Morph adjusted for the influence of depth.
Color Morph and Dominant Symbiont per Bay Area
Latitude=aggregate(Latitude~Bay.Area, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Bay.Area, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Bay.Area", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
merged$ColDom <- interaction(merged$Color.Morph, merged$Dom)
model <- multinom(ColDom~newDepth+Bay.Area, merged)
means <- lsmeans(model, specs = "Bay.Area")
means
newdat <- data.frame(Bay.Area = levels(merged$Bay.Area),
newDepth = mean(merged$newDepth, na.rm=T))
pred <- predict(model, newdata = newdat, "probs")
new <- data.frame(Bay.Area=as.character(newdat[,1]), pred)
XY4 <- merge(XY, new, by="Bay.Area", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Bay.Area
par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)
XY4$Bay.Area <- as.numeric(XY4$Bay.Area)
apply(XY4, MARGIN=1, FUN=function(reef) {
floating.pie(xpos = reef["X"], ypos = reef["Y"],
x=c(reef["Brown.C"], reef["Orange.C"], reef["Brown.D"],
reef["Orange.D"]), radius = 20,
col = c("blue", "cornflowerblue", "red", "lightpink"))
})
legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("blue", "cornflowerblue", "red", "lightpink"))
par(new=T, mar=c(14,21,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.87, 21.41, -157.75, 21.52)
box()
par(mar=c(3,3,0,0))
plot(HI, xlim=c(-157.82, -157.82), ylim=c(21.42, 21.51),
lwd=2, col="gray")
legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"),
fill = c("turquoise4", "steelblue1", "tomato", "coral"), bg = "white")
add.pie(z = c(42, 23, 5, 30), x = -157.8109, y = 21.46222, labels = NA, radius = 0.005,
col = c("turquoise4", "steelblue1", "tomato", "coral"))
add.pie(z = c(36, 33, 2, 29), x = -157.8286, y = 21.47748, labels = NA, radius = 0.005,
col= c("turquoise4", "steelblue1", "tomato", "coral"))
add.pie(z = c(40, 24, 8, 28), x = -157.7965, y = 21.44077, labels = NA, radius = 0.005,
col= c("turquoise4", "steelblue1", "tomato", "coral"))
box()
par(new=T, mar=c(17,24,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83")
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83"))
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.87, 21.41, -157.75, 21.52)
box()